using NUnit.Framework;
using dnAnalytics.LinearAlgebra;
using dnAnalytics.LinearAlgebra.Solvers.Preconditioners;

namespace dnAnalytics.UnitTests.LinearAlgebra.Solvers.Preconditioners
{
    [TestFixture]
    [Category("Managed")]
    public sealed class IlutpPreconditionerTest : PreconditionerTest
    {
        private double mDropTolerance = 0.1;
        private double mFillLevel = 1.0;
        private double mPivotTolerance = 1.0;

        [SetUp]
        public void Setup()
        {
            mDropTolerance = 0.1;
            mFillLevel = 1.0;
            mPivotTolerance = 1.0;
        }

        private SparseMatrix CreateReverseUnitMatrix(int size)
        {
            SparseMatrix matrix = (SparseMatrix)MatrixBuilder.CreateMatrix(size, MatrixType.Sparse);
            for (int i = 0; i < size; i++)
            {
                matrix[i, size - 1 - i] = 2;
            }
            return matrix;
        }

        private IlutpPreconditioner InternalCreatePreconditioner(SparseMatrix matrix)
        {
            IlutpPreconditioner result = new IlutpPreconditioner(matrix);
            result.DropTolerance = mDropTolerance;
            result.FillLevel = mFillLevel;
            result.PivotTolerance = mPivotTolerance;
            return result;
        }

        internal override IPreconditioner CreatePreconditioner(SparseMatrix matrix)
        {
            mPivotTolerance = 0;
            mDropTolerance = 0.0;
            mFillLevel = 100;
            return InternalCreatePreconditioner(matrix);
        }

        protected override void CheckResult(IPreconditioner preconditioner, Vector vector, Vector result)
        {
            Assert.AreEqual(typeof(IlutpPreconditioner), preconditioner.GetType(), "#01");
            IlutpPreconditioner factorization = preconditioner as IlutpPreconditioner;
            // Get the Upper triangular matrix U
            // Get the Lower triangular matrix L
            // Get the diagonal matrix D
            // M = (L * U)
            // Compute M * result = product
            // compare vector and product. Should be equal
            Matrix l = factorization.LowerTriagonalMatrix();
            Matrix u = factorization.UpperTriagonalMatrix();
            // Get the pivot matrix
            int[] pivots = factorization.Pivots();
            SparseMatrix p = (SparseMatrix)MatrixBuilder.CreateMatrix(l.Rows, MatrixType.Sparse);
            for (int i = 0; i < p.Rows; i++)
            {
                p[i, pivots[i]] = 1.0;
            }
            Matrix temp = l.Multiply(u);
            Matrix matrix = temp.Multiply(p);
            Vector product = VectorBuilder.CreateVector(result.Count, VectorType.Dense);
            matrix.Multiply(result, product);
            for (int i = 0; i < product.Count; i++)
            {
                Assert.AreEqual(vector[i], product[i], epsilon, "#02-" + i.ToString());
            }
        }

        [Test]
        public void SolveReturningOldVectorWithoutPivoting()
        {
            int size = 10;
            SparseMatrix newMatrix = CreateUnitMatrix(size);
            Vector vector = CreateStandardBcVector(size);
            // set the pivot tolerance to zero so we don't pivot
            mPivotTolerance = 0.0;
            mDropTolerance = 0.0;
            mFillLevel = 100;
            IPreconditioner preconditioner = CreatePreconditioner(newMatrix);
            preconditioner.CreatePreconditioner();
            Vector result = VectorBuilder.CreateVector(vector.Count, VectorType.Dense);
            preconditioner.Solve(vector, result);
            CheckResult(preconditioner, vector, result);
        }

        [Test]
        public void SolveReturningOldVectorWithPivoting()
        {
            int size = 10;
            SparseMatrix newMatrix = CreateUnitMatrix(size);
            Vector vector = CreateStandardBcVector(size);
            // Set the pivot tolerance to 1 so we always pivot (if necessary)
            mPivotTolerance = 1.0;
            mDropTolerance = 0.0;
            mFillLevel = 100;
            IPreconditioner preconditioner = CreatePreconditioner(newMatrix);
            preconditioner.CreatePreconditioner();
            Vector result = VectorBuilder.CreateVector(vector.Count, VectorType.Dense);
            preconditioner.Solve(vector, result);
            CheckResult(preconditioner, vector, result);
        }

        [Test]
        public void CompareWithOriginalDenseMatrixWithoutPivoting()
        {
            SparseMatrix sparseMatrix = (SparseMatrix)MatrixBuilder.CreateMatrix(3, MatrixType.Sparse);
            sparseMatrix[0, 0] = -1;
            sparseMatrix[0, 1] = 5;
            sparseMatrix[0, 2] = 6;
            sparseMatrix[1, 0] = 3;
            sparseMatrix[1, 1] = -6;
            sparseMatrix[1, 2] = 1;
            sparseMatrix[2, 0] = 6;
            sparseMatrix[2, 1] = 8;
            sparseMatrix[2, 2] = 9;
            IlutpPreconditioner ilu = new IlutpPreconditioner(sparseMatrix);
            ilu.PivotTolerance = 0.0;
            ilu.DropTolerance = 0;
            ilu.FillLevel = 10;
            ilu.CreatePreconditioner();
            SparseMatrix l = ilu.LowerTriagonalMatrix() as SparseMatrix;
            // Assert l is lower triagonal
            for (int i = 0; i < l.Rows; i++)
            {
                for (int j = i + 1; j < l.Rows; j++)
                {
                    Assert.AreEqual(0.0, l[i, j], epsilon, "#01-" + i.ToString() + "-" + j.ToString());
                }
            }
            SparseMatrix u = ilu.UpperTriagonalMatrix() as SparseMatrix;
            // Assert u is upper triagonal
            for (int i = 0; i < u.Rows; i++)
            {
                for (int j = 0; j < i; j++)
                {
                    Assert.AreEqual(0.0, u[i, j], epsilon, "#02-" + i.ToString() + "-" + j.ToString());
                }
            }
            SparseMatrix original = l.Multiply(u) as SparseMatrix;
            for (int i = 0; i < sparseMatrix.Rows; i++)
            {
                for (int j = 0; j < sparseMatrix.Columns; j++)
                {
                    Assert.AreEqual(sparseMatrix[i, j], original[i, j], epsilon, "#03-" + i.ToString() + "-" + j.ToString());
                }
            }
        }

        [Test]
        public void CompareWithOriginalDenseMatrixWithPivoting()
        {
            SparseMatrix sparseMatrix = (SparseMatrix)MatrixBuilder.CreateMatrix(3, MatrixType.Sparse);
            sparseMatrix[0, 0] = -1;
            sparseMatrix[0, 1] = 5;
            sparseMatrix[0, 2] = 6;
            sparseMatrix[1, 0] = 3;
            sparseMatrix[1, 1] = -6;
            sparseMatrix[1, 2] = 1;
            sparseMatrix[2, 0] = 6;
            sparseMatrix[2, 1] = 8;
            sparseMatrix[2, 2] = 9;
            IlutpPreconditioner ilu = new IlutpPreconditioner(sparseMatrix);
            ilu.PivotTolerance = 1.0;
            ilu.DropTolerance = 0;
            ilu.FillLevel = 10;
            ilu.CreatePreconditioner();
            SparseMatrix l = ilu.LowerTriagonalMatrix() as SparseMatrix;
            SparseMatrix u = ilu.UpperTriagonalMatrix() as SparseMatrix;
            int[] pivots = ilu.Pivots();
            SparseMatrix p = new SparseMatrix(l.Rows);
            for (int i = 0; i < p.Rows; i++)
            {
                p[i, pivots[i]] = 1.0;
            }
            SparseMatrix temp = l.Multiply(u) as SparseMatrix;
            SparseMatrix original = temp.Multiply(p) as SparseMatrix;
            for (int i = 0; i < sparseMatrix.Rows; i++)
            {
                for (int j = 0; j < sparseMatrix.Columns; j++)
                {
                    Assert.AreEqual(sparseMatrix[i, j], original[i, j], epsilon, "#01-" + i.ToString() + "-" + j.ToString());
                }
            }
        }

        [Test]
        public void SolveWithPivoting()
        {
            int size = 10;
            SparseMatrix newMatrix = CreateReverseUnitMatrix(size);
            Vector vector = CreateStandardBcVector(size);
            IlutpPreconditioner preconditioner = new IlutpPreconditioner(newMatrix);
            preconditioner.PivotTolerance = 1.0;
            preconditioner.DropTolerance = 0;
            preconditioner.FillLevel = 10;
            preconditioner.CreatePreconditioner();
            Vector result = VectorBuilder.CreateVector(vector.Count, VectorType.Dense);
            preconditioner.Solve(vector, result);
            CheckResult(preconditioner, vector, result);
        }
    }
}